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Abstract 


We formulate the Dirac equation for a massive neutral spin-half particle on a rotating black 
hole spacetime, and we consider its (quasi)bound states: gravitationally-trapped modes which are 
regular across the future event horizon. These bound states decay with time, due to the absence of 
superradiance in the (single-particle) Dirac field. We introduce a practical method for computing 
the spectrum of energy levels and decay rates, and we compare our numerical results with known 
asymptotic results in the small-M/r and large-Mp regimes. By applying perturbation theory in 
a horizon-penetrating coordinate system, we compute the ‘fine structure’ of the energy spectrum, 
and show good agreement with numerical results. We obtain data for a hyperfine splitting due to 
black hole rotation. We evolve generic initial data in the time domain, and show how Dirac bound 
states appear as spectral lines in the power spectra. In the rapidly-rotating regime, we find that 
the decay of low-frequency co-rotating modes is suppressed in the (bosonic) superradiant regime. 
We conclude with a discussion of physical implications and avenues for further work. 
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I. INTRODUCTION 


Einstein’s theory of general relativity (1915) and Dirac’s relativistic wave eqnation (1928) 
are cornerstones of physics and ‘two households alike in dignity’. Although Dirac may have 
remarked that ‘it is more important to have beauty in one’s equations than to have them 
fit experiment’, these theories have also endured ever-more-stringent experimental tests. 

Thus far, there has been little opportunity to test the beautiful theory of fermions on 
curved spacetimes, developed by Schrodinger, Bargmann, Pauli [l| (1933) and others [3, T 


except for in a restrictive ‘Newtonian’ setting |4,l5|. This may change in forthcoming decades. 
Black holes - a radical consequence of Einstein’s theory which appear to play a pivotal role 
in the development of galactic structure - may act as crucibles for stringent tests of strong- 
held phenomena j^, and the interaction of fermionic helds with black holes, the topic of this 
article, is of interest from a range of perspectives. 

In standard experience, there is a vast difference between the typical gravitational and 
quantum-mechanical length scales. A gravitating mass M interacting with a quantum held 
of mass /i is characterized by a dimensionless parameter 

Mp, GMfj, Th 
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where is the horizon radius of the gravitating mass, Ac is the Compton wavelength of the 
held, and mpi is the Planck mass. As is typically measured in kilometres, and Ac in pico 
or femtometres, it is clear that Mp ^ 1 in standard astrophysical scenarios (henceforth we 
will typically omit dimensionful constants, setting G = h = c = mp\ = 1). Yet one may 
also envisage scenarios in which this is not the case, e.g., for (i) standard-model helds bound 
to ‘light’ primordial black holes, and (ii) ultra-IDht helds (arising from e.g. string theory 
compactiheations) bound to galactic black holes 7|, l8| . It has been suggested that ultra-light 
helds could provide a resolution of the dark matter problem (see e.g. Refs. ilul). 


The ‘no-hair’ (Israel-Carter) conjecture 12 h 14| suggests that, once perturbed, a black hole 


will rapidly return to a stationary state, changing only a small number of physical parameters 
in the process: its mass M, angular momentum J, and charge Q. The conjecture has 
been elevated to a theorem for various scenarios involving massless scalar, electromagnetic 


and gravitational helds 12l4l7l|. However, if a perturbing held is endowed with mass, the 


quantum/classical correspondence principle suggests that the ‘no-hair’ picture is incomplete. 
Outside the innermost stable circular orbit {visco = QGM/c^ for Schwarzschild black holes), 
a compact body of mass p may orbit a black hole indehnitely, at least up to a gravitational 
radiation-reaction timescale Trad ~ (M/p)(GM/c^). In the semi-classical regime Ac <C rh 
(i.e. Mp 3> 1) a straightforward WKB analysis (see Appendix |A]) suggests that a massive 
held possesses bound states E < pc^ localized around orbiting timelike geodesics. For 
circular orbits in the weak-held regime {tq S> M), 
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where the first term is the rest mass energy, the second term is the ‘Newtonian’ binding 
energy, and the third term gives the energy-level spacing in terms of the orbital frequency 
ffo = ^/GM/ry^. In cases where KIq -C jjLC? the discrete spectrum becomes indistinguish¬ 
able from a continuum. 

There is a subtlety in the above picture. Localized states are (generically) asymmetric 
under time-reversal, due to non-Hermiticity in the Dirac equation on black hole spacetimes. 
Flux is absorbed through the event horizon, which acts as a one-way membrane. Notwith¬ 
standing, lifetimes become exponentially-long when an angular momentum barrier separates 
the orbit from the event horizon. That is, we have r ~ where /3 is some tunnelling 

coefficient which depends only on the ratio of Mp,/(£ -1- 1/2), with £ the orbital angular 
momentum number. 

How does a black hole’s rotation affect this picture? Consider a Kerr black hole, with 
angular momentum aM. It is well-established that low-frequency modes of bosonic helds 
may be amplified by a process known as superradiance. Localized bosonic states within 
the superradiant regime ojoj < 0, where oj = oj — mfln with Qh = a/{2Mrh) the angular 
frequency of the horizon, and m the azimuthal number of the mode, may grow with time, 
cau sing a superradiant instability to develop: this is the so-called ‘black hole bomb’ scenario 
18|,ll9j. Although growth, like decay, is exponentially-suppressed in the semi-classical regime 
Mfi S> 1 (with e-folding time r > 10'^exp(1.84Mp) 20|), it becomes significant for M/i ~ 
0(1) (with Tmin ~ 5.88 X lO^GM/c^ for the i = m = 1 mode of the scalar held at Mp = 0.45 
and a = 0.997M 2l424|). By contrast, the (single-particle) Dirac equation is not subject 
to superradiance 25l-l29j| . and thus all states decay. See Ref. 30[ for a recent summary of 
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for recent work on the evolution 


superradiance and associated phenomena, and Refs, 
of superradiant instabilities. 

Many authors [ 3 - 0 , 0 - 0 ] have considered the ‘bound states’ of massive helds on 
Schwarzschild, Kerr and Kerr-Newman spacetimes. A range of terminology has been used, 
including ‘quasiresonances’, ‘quasilevels’ |52L ‘quasistationary states’, ‘quasibound states’ 
j3 |. ‘dynamical resonance states’ 53], ‘wigs’ 54] and ‘graviatoms’ 48]; here we use ‘bound 
states’ generically, without any implication of time-reversal symmetry. Bound states are 
not to be confused with ‘quasinormal modes’, which are radiative in character. Loosely 
speaking, quasinormal modes are associated with a maximum (rather than a minimum) in 


the ehective potential; they were studied in the massless Dirac case in Refs. [SSHS?]. 


Regardless of diherences in terminology, in the regime Mp/(£-1-1/2) < 1, a unihed picture 
of bound states emerges. Let oj = ujr + iui denote the real and imaginary parts of the mode 
frequency, with uj < 0 {uj > 0) corresponding to decay (growth). To lowest order, the 
spectrum of spin-zero and spin-half particles is hydrogenic, with 

(Mp)2 


^ 1 


2n2 


(3) 


with n G 1,2,... the principle quantum number. In this regime, the bound states have a 
typical radius of noo, where Oq = {Mfi)~‘^{GM/c^) is the gravitational ‘Bohr radius’. The 


3 





























imaginary part of a bosonic field of spin s (scalar s = 0, Proca s = 1 and massive graviton 


s = 2) scales according to [34 
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where S G {—s, —s -|- 1,..., s} is the spin projection. By comparison, for the Dirac held it 
was shown in Ref. [37| that, for a M, 


cu/Z/i oc —(M/r) 


4Z+2S+4 


( 5 ) 


where S G { — 1/2,1/2} is the spin projection. This highlights a key point: in the regime 
/i < mflH, bosonic modes become unstable, due to superradiant growth, whereas fermionic 


modes remain stable 
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A more subtle point here is that the index in Eq. ([5]) is smaller 


by 1 than would be anticipated from Eq. (|1]). 

In this paper we formulate a practical method for computing Dirac bound states on the 
Kerr spacetime. Our aim is to use highly-accurate numerical data to verify and test various 
asymptotic results, such as those above; and to investigate the rich phenomenonology of the 
intermediate regime Mp 0(1). 

Rather remarkably, the massive Dirac equation on the Kerr spacetime admits a complete 
separation of variables. This reduces our problem to the analysis of (coupled) ordinary 


differential equations. Separability was hrst shown by Unruh 62| in 1973 for the massless 


held, and by Chandasekhar 63|, in 1976 for the massive held. Whereas Unruh’s analysis 
used the standard 4-spinor formalism, Chandrasekhar’s analysis employed the Newman- 
Penrose 2-spinor formalism. In 1979, Carter & McLenaghan 65| showed that the Dirac 
operator commutes with a ‘generalized total angular momentum operator’ constructed from 
the Killing-Yano tensor. In 1984, Kamran & McLenaghan 66| described (see theorem 3) a 


class of spacetimes for which the Dirac equation admits a separation of variables in the Weyl 
representation. In 1993, McKellar et al. 67| conducted an explicit separation of variables 


in t 


le 4-spinor formalism, using the Pauli-Dirac representation. Finster and collaborators 


68l-l7l|. and others 72|, l73| have also used the 4-spinor formalism, along with the Weyl 


representation. 

From one point of view, superradiance is a necessary consequence of the second law 
of black hole thermodynamics: the horizon area A is a non-decreasing function of time. 
How, then, is the absence of superradiance consistent with this law? A key assumption 
underpinning the second law is that the weak energy condition holds, i.e., > 0 for 

any timelike vector field, < 0, where T„„ is the stress-energy tensor. This condition is 


violated by the (single-particle) Dirac held 621. 

This paper is organised as follows. In Sec. Hi] we review the formalism for describing a 
neutral spin-half held on a Kerr black hole spacetime. In Sec. IIIII we describe our methods 
for calculating the bound state spectrum. In Sec. jTVjwe explore the bound state spectrum on 
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Schwarzschild and Kerr spacetimes, developing and testing asymptotic resnlts in the regimes 
Mjj, -C 1 and M/r 1. In Sec. |V]we discnss onr findings. 

Thronghont, we set G = c = h = 1 and adopt a metric signatnre +2 (except for in 
Appendix |D]). Coordinate indices are denoted with Greek letters a, /5 ,7 ... and tetrad basis 
indices with Roman letters a,b,c,.... Symmetrization and antisymmetrization of indices is 
denoted with ronnd and sqnare brackets, () and [], respectively. We nse 9^, and to 
denote partial, covariant and spinor derivatives, respectively. 


II. FORMALISM 

In this section we ontline the formalism for describing spin-half fields on the Kerr black 
hole spacetime (N.B. expert readers may wish to proceed immediately to Sec. Illlj) . 


A. Kerr spacetime 


1. Coordinate systems 


The region ontside the event horizon of a Kerr black hole may be described with the 
Boyer-Lindqnist coordinate system {t,r,9,(j)}, in which the line element ds^ = g^^dx^dx^ 
takes the form 


ds^ 


2Mr\ , 9 
1 -^ dt^ 




V 


AaMr sin^ 6 ^ , , , , 9 

-dtdct) + ^dr^ 




A 


2 1 2 2 2 Mra^sin^P\ , o „ , ,0 

+p^d9^ H-^- ) sm^ 


( 6 ) 


where + a? cos^ 9 and A = — 2Mr -|- a?. The Kerr spacetime has two horizons 

at = M ± -sj— a^, and two stationary limit snrfaces, at rs± = M ± -v/M^ — cos^ 9. 
The angnlar velocity of the event horizon is Qh = a/2Mr^. 

A deficiency of the Boyer-Lindqnist system is that it takes an infinite coordinate time 
t for ingoing geodesics to reach the onter horizon at r = r+. To properly describe snch 
geodesics, we may employ the ‘ingoing Kerr’ coordinate system {t, f, 9, 0} with r = r, 9 = 9 

d^ = d(j) + 


ana 


~ 2Mr 

dt = dt -\ --— dr, 

A 


Explicitly, t = t + a{r), cj) = (j) + (3{r) where 


a(r) 


0 (r) 


2Mr+ 

-In r — r+ 
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a , r — r+ 

-In - 
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In this coordinate system the inverse metric takes a simple form, 


/-(p2 + 2Mr) 2Mr 0 0 \ 
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( 10 ) 


and the ingoing principal null geodesics are simply given by dr = —dt. In either coordinate 
system, the metric determinant is given by \/—g = sin 9. 


2. Canonical tetrad 


Let us introduce a tetrad of vectors e" = {cq , e", , e^} which we take to be an orthonor¬ 

mal basis, i.e., gape^e^ = rjajs where rjab = diag(—1,1,1,1). Roman and Greek letters are 
used for tetrad and coordinate indices, respectively. Roman indices are raised and lowered 
with pab- It follows that e“ = g^'^gapeb gap = VabeaGp. 

We will employ the ‘canonical’ orthonormal (symmetric) tetrad for the Kerr spacetime 


introduced by Carter [74, 


75 


VIZ., 


-I- 


eo = 


P\/A 

: = ^/p. 

asin0 


^0 ~ 


ei = 


P 




= i/p, 


Gn = 


p sin 9 ’ 


( 11 ) 


with inverse components given in flBlD such that ds^ = pab{Glj_dx^){eldx'^). It is straightfor¬ 
ward to find the components of the canonical tetrad in the ingoing Kerr coordinate system. 


using g; = 


3. Spacetime symmetries and conservation laws 

The spacetime admits two independent Killing vectors, = {dt)^ and with 

the defining properties = 0 = Furthermore, the spacetime admits a Killing-Yano 

tensor with the defining properties = fi^^] and namely 

= 2 (a cos 9 + r . (12) 

Now consider the stress-energy tensor associated with the field, which is symmetric 
in its indices and divergence-free = 0). With the Killing vectors, we may form 

two independent divergence-free vectors, and 

associated with energy and azimuthal angular momentum, respectively. A third divergence- 
free vector is given by the Dirac probability current (see Secs. Ill B 31 and IIIBSp . 
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From each divergence-free vector (V^J^ = 0) one may obtain a conservation law via 
Gauss’ theorem, 






'dV 


111 


(13) 


where V is a contiguous four-volume bounded by a three-volume (hypersurface) dV. Here, 
the volume element is dV = \/—g dtdrdOdcp, and the 3-surface element is dehned in 
terms of the metric induced on the boundary hypersurface in the standard way 76|. We 
may construct a four-volume of inhnitessimal extent At conhned between twin spacelike 
hypersurfaces at f = to and t = to + At and twin timelike hypersurfaces r = ri and 
r = r 2 (where to and r 2 > ri > r+ are constants). This construction leads to a quasi-local 
conservation law in the form 

d 


dt 


'ri 


{p^r)dndr )■ = - 


r2 


{p^r)dn 


(14) 


ri 


where dG = sin OdOdcj). Similarly, by considering hypersurfaces of constant t, we may obtain 
the corresponding expression in the ingoing-Kerr coordinate system. 


di 



(pV')dGdr 



{p^r)dn 


r2 

ri 


(15) 


where J* = J* -|- J'" and = J'’. Note that the t = const and i = const hypersurfaces 

are rather different in character, as the former approaches the bifurcation point, whereas the 
latter penetrates the future horizon. Physically, we expect (and hnd) that the ingoing-Kerr 
quantity J* is hnite as r —>■ r_|_, whereas the Boyer-Lindquist version J* is not. 


B. Spin-half fields 

1. The Dirac equation 

The Dirac equation on a curved spacetime (in signature —h -|--|-) takes the form 

- pj T = 0, (16) 

where p is the held mass, T is a Dirac four-spinor, are Dirac four-matrices, and the spinor 
covariant derivative Dy is 

D, = d,- T,. (17) 

The spinor connection matrices are dehned, up to an additive multiple of the unit matrix, 
by the relation 

- r.7^ + 7'^r, = o. (is) 

where is the affine connection. A suitable choice satisfying (flSj) makes use of the spin 
connection Uahc (described below), 

Ba = (19) 
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Here, 7 “ and 7 “ denote sets of 4 x 4 matrices satisfying the anticommutation relations 

{7“, 7^} = {r, 7'} = 27“'/4, (20) 

where {A, B} = AB + BA. Note that the former set 7 “ are functions of spacetime position, 
whereas the latter 7 “ have constant components. The two sets may be related with any 
orthonormal tetrad, 

7“ = • (21) 

We will employ the canonical tetrad of Eq. (nj. Matrix indices are raised/lowered in the 
standard way: ja = VabJ^ and 7 q, = ■ 


2. The spin connection 


The spin-connection u^ab arises naturally when one considers how a generalized covariant 
derivative should act upon tensors with mixed coordinate and basis components. We 

start with 

= d,Al + - T\^A:. (22) 

Imposing metric compatibility londs to the usual definition for the affine 

connection F" . Imposing tetrad compatibility — 0) leads to 

‘r'/t = - 4d,et. (23) 

To obtain the spin connection without hrst calculating the affine connection one may use 

^fiab 2^^ (.^abc T ^cab ^bca) t (24) 

where Xabc = e.^{dve-bfi — dfj_eb^) e'^ (see Ref. for details). The spin connection for the 
canonical tetrad is listed in Appendix |Bl 


3. Current and stress-energy 


The Dirac current is given by 

J^^ = (25) 

where 4/ = 4/la is the Dirac conjugate, with denoting the usual Hermitian conjugate. 
The Hermitizing matrix a must satisfy the conditions 


ay 


+ 7 ^^q; = 0, d^a + Tj^a -|- aV^ = 0. 


We choose a = — 7 ° [77 1. 

The (symmetric) stress-energy is given by 


T = - 


T7(^D,)T - 7,)T 


(26) 


(27) 


where the covariant derivative of the conjugate spinor is = 9^4/ -|- 4/r^. The Dirac 
current and stress-energy are covariantly conserved, = 0 = 





Matrix representation 


We will use the Weyl/chiral representation, in which 


7° = 


O I 
I o 


7 = 


O ai 
-CTi O 


i = 1,2,3, 


(28) 


where I is the 2x2 identity, O is the 2x2 zero matrix, and are the Pauli matrices 


0 1 


0 -i 


' 1 0 / ’ \ z 0 


= 


1 0 
0 -1 / ’ 


(29) 


which satisfy |{cTj, aj} = 6ijl + ieijkUk- 

In fact, we will choose our matrices 7 * (Sec. Ill B ip to be a cyclic permutation of the set 
7 * multiplied by the unit imaginary, as follows: 7 ^ = Z 7 ^, 7 ^ = * 7 ^, 7 ^ = * 7 ^ and 7 ° = * 7 °. 
Note that making the permutation is equivalent to relabelling the tetrad legs; alternatively 
one may take the view that we are applying a (constant) similarity transform 7 ^ —)■ S'y^S~^, 
'h —>■ S'd'. Including the factor of i is necessary in order to satisfy anticommutation relations 
(| 20 l) on a spacetime of positive signature. 

V-' 


We will write the Dirac four-spinor as d' = 


>+. 


where and 'll)- are (left- and right- 


handed) two-spinors, which may be projected out from d' with the operators ^ (J ± 75 ) 
where 

o' 

O I 


7 ^ = Z7°7^7^7^ = 


(30) 


5. Separation of variables 


Let us introduce the complex quantity g = r + ia cos 6, and its conjugate g* = r — ia cos 6, 
such that = gg*. A short calculation (using Eq. (IT^ and Appendix |B]) shows that the 
spin-connection matrices P^ take the form given in Eq. flB3p . Next we obtain 




1 I O s^ai s*a3 

2 [-(s^ai + SrO-a) O 


(31) 


where Sr = and sg = e\ -Y^-§q {gsinO). This result for the spin-connection 

matrix suggests a natural ansatz for the wavefunction. 




^-1/4 


g 

£>- 7 +) 


(32) 
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where ri± are two-spinors. Multiplying the Dirac equation fflB]) by 
yields a pair of two-spinor equations, 

1 


( 


.1/2 Q 


v ^ 


n /2 


Va 


[(r^ -F a^)dt ad^] + a^\fAdr + ai {de \ cot 6) + (a sin 6dt + esc 6d^) > p± 


^ip (r = 1 = ia cos 6) r/=p = 0. 

(33) 


Now, separating harmonic temporal and azimuthal dependence with 

where m is a half-integer, leads to a remarkable separation of the r and 9 parts, 

+ am] + (T 3 \/Ac?,.| =F iprp± 

+ {<Ji [de + ^ cot 9) + 102 (—aoj sm9 + m esc 6*)} i/=p -|- ap cos 9p± =0, 
(cf. Eq. (8) in Ref. j^; see also Appendix [C]) . Introducing the ansatz 


(34) 


(35) 


V+ = 


Ri{r)S,{9) 

Mr)S2{9))^ 


V- = 


R2{r)S,{9) 
Ri(r)^2(0)/’ 


(36) 


and multiplying by as leads to twin pairs of coupled first-order ordinary differential equa¬ 
tions. 


\/A {dr — iK/A) Ri = {\ + ipr) R 2 , 
\/A {dr + iK/A) R 2 = {X — ipr) Ri, 


and 


de + - cot 9 — m esc 6* -f ao; sin 6* ) Ri = (-I-A -|- ap cos 9) S' 2 , 


de + - cot 9 + mcsc9 — aujsin9 ) S '2 = (—A -|- ap cos 9) Si, 


(37) 

(38) 


(39) 

(40) 


where K = {r^ + aS)uj — am and A is the separation constant. Eqs. fl37)l - fl40|) are equivalent 
to the coupled ordinary differential equations originally obtained by Chandrasekhar 63. Im|. 


6. Angular solutions 

In the non-rotating case (a = 0), the eigenvalues of the angular equations fl3^ -fl40 |) are 
integers A G Z\0 = {..., —2, —1, -|-1, 2 ...}, and the solutions are spin-weighted spherical 
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harmonics (see Sec. 3 in Ref. 78|). The total angular momentum, j, a half-integer, is related 
to the a = 0 eigenvalue in a simple way: A = 'P{j -|- 1/2), where V = ±1. We may dehne 
the orbital angular momentum, £, an integer, via = j + \V. 

In the rotating case the solutions of Eqs. fl3^ -fH0 |) are known as mass-dependent 
spheroidal harmonics (see Ref. [78| for a review). We let Si = _S'a and S 2 = +S\, 
where A = {j, m, V, au, a/i} and we note the key symmetry relation 


±Sji{9)=V{-iy+^^SA{n-9), 

and normalization condition 



(41) 


(42) 


7. Radial solutions 


Let R = I . In the near-horizon region r r+, the ‘horizon-ingoing’ radial solution 
R 2 


has the asymptotic form 


where 


Rhor ~ I ^ I exp 


U! = u — mflff, 

A -|- i/ir+ 




r_|_ — r_ — 4zMa)r_|_ 

Here r* is the tortoise coordinate dehned by 

dr* _ (r^ -F a^) 
dr A 

The ‘horizon-outgoing’ solution is obtained by the interchange Ri —)• R 2 and R 2 ^ R\- 
In the far-held region, r —)■ cxo, the propagating solutions take the form 


(43) 

(44) 

(45) 

(46) 


R) 


^±7g±ipr 

A/2p(ar^^ 


-A 


+ 0{r ) 


7 = 


- /i 2 ) 

P 


(47) 


where p = — p?. The bound solutions {u < p) may be obtained by replacing ip with 

q = — in the above. 

Multiplying (ITni and fl38l) by Rl and R 2 , respectively, and taking the difference, leads to 
the Wronskian relationship A _ IR 2 P) = 0. Now let us consider a modal solution of 

the form 


AfeRhor AqqRqq -I- Aq^Rqq, 


(48) 
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where Ah and are complex constants. From the constancy of the Wronskian, it follows 
immediately that 

(49) 

that is, 7^ = 1 — T, where the reflection and transmission coefficients are TZ = 1^4+/ \A'^\^ 
and T = \Ah!\^ / As T > 0, it is clear that superradiance does not occur in the modal 

reflection/transmission problem. 


8. The Dirac current and absence of superradiance 


We now insert fl36|l (see also fl3^ and (El])) into the dehnition of the Dirac current 
The radial component is given by 

^ (l^d - l-R'it) (|sy + |S2|"). 


(60) 


Let us refer back to the conservation law fflTD . The integral of the radial current evaluated 
at the lower limit ri = r_|_ has a natural interpretation as the flux passing into the horizon; 
we have 

where N is the number density. We saw in the previous section that the right-hand side is 
negative for horizon-ingoing solutions: hence dN/dt < 0. This provides further confirmation 


that superradiance is absent [25 


27|,l62|. 


Next, we obtain expressions for the components of the Dirac current in the canonical 
basis. 


J(0) ^ gO J/4 _ 

^ P\/A 

j(3) ^ gS jM = J_ 

^ p\/A 


Alm{RlR2)Re{SlS2), 


(52) 

(53) 


The temporal component of the current, .P = is 


P = 


-I- P 


+ \R2?) (|5i|2 + 152^ + ^^^5^Im(i?*i?2) Re(5*52) 


(54) 


Note that the factor of A in the denominator implies that P diverges in the limit r —)■ r+. 
This is due to the coordinate singularity in the Boyer-Lindquist system. We may instead 
consider the temporal component in the ingoing-Kerr coordinate system, given by P = 

-(E])), i.e. 


P -I- Ji" evaluated at t = t — a(r), 0 = 0 — /3(r) (cf. Eq. 


P = 




IR 2 P + (r^ + 2Mr + P) lA-^^^ipj (|5ip + |^ 2 p) 
-|-4a sin 9 lm{A~^f‘^R\R 2 ) Re(S'J' S 2 ) 


( 55 ) 
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where 




fc = 1,2. 


(56) 


From the asymptotic form of the ingoing solution Rhor, Eq- H3|) . it is clear that A ~ 

0 ( 1 ) and thus is finite on the (future) horizon r = r+, as expected. 


9. Violation of weak energy principle 


The weak energy condition states that > 0 for any timelike vector t^. Let 

us now introduce the quantity S = —T^^e^CQ, noting that Cq is indeed timelike. A short 
calculation gives S = Sq + Si, where 


50 = [(r^ + a^)cu - am] (|Rip + jRaP) (|5'ip + 

51 = 5^ [cos0 (|Ri|2 + IR 2 P) (l^ir - 1 *^ 21 ') 

Ip 

4r 

- sin e-= Re{RlR 2 ) Re(^*^ 2 ) • 

V A 


(57) 


(58) 


The hrst term Sq arises from the partial derivatives in Eq. fE7)) . and the second term Si arises 
from the spin-connect ion matrices. The hrst term dominates over the second as A —)■ 0. It 
is clear that Sq is not positive-dehnite; it is negative in the near-horizon region if ca < mfln- 
Thus the weak-energy condition is violated for the Dirac equation on Kerr spacetime in 
the superradiant regime, at least at the ‘one particle’ level (for second quantization, see 

Ref. 0 )- 


III. METHODS 


In this section we outline practical frequency-domain (Sec. IIII Ap and time-domain 
(Sec. IIIIBp methods for investigating the bound state spectrum. 


A. Frequency-domain method 


We now formulate a practical method for computing the discrete spectrum of bound 
states. In essence, we wish to solve the radial equations (1Tni - (]38|l subject to certain boundary 
conditions. On the horizon, the solution should be hnite in the Kerr-ingoing coordinate 
system; towards inhnity, the solution should decay exponentially. Thus, the integral of the 
probability density fl5^ should be hnite across the exterior region r > r_|_. 

A powerful method for computing the spectrum of quasi-normal modes was introduced in 


Ref. 


, in the context of massless helds. With an ansatz adapted to the boundary condi¬ 


tions, the differential equation for massless bosonic helds generates a three-term recurrence 
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relation. The minimal solution of the recurrence relation is found by solving a continued- 
fraction equation. This approach was adapted to compute the bound state spectrum in 
Ref. j^. We follow this approach below. 


1. Angular solutions 


Three-term recurrence relations for the angular equations 


tained in Refs. 78 


81 


were previously ob- 
82| . We follow the method of Ref. TSj], in which the mass-dependent 


spin-half spheroidal harmonics are decomposed in a basis of spin-half spherical harmonics, 
leading to a three-term recurrence relation for the expansion coefficients {bk} (where fc is a 
half-integer): 


c^kbk+i + Pkbk — 0 , 
akbk+i + jSkbk + Akbk-i = 0 , 


k = \m\, 

k = \m\ A-1 , \m\ A 2, 


where 


CXk — (op- + 


\/{k A 1)^ — m? 


l3k = ek{kAl/2)[l- 


2(fc + l) 
amoj 


A 


ajum 


k{k Al) J 2k{k A 1) 


-A, 


7 fc = (ap - ekau) 


yjk"^ — -m? 
2fc 


(59) 

(60) 


(61) 

(62) 

(63) 


and Cfc = with j = i-A V12. The boundary conditions at the poles are satished 

by the minimal solution of the recurrence relation, which is found by solving a continued- 
fraction equation for A, 


Q Q^|m|/3|m|+l Q^|m.|+ld|m|-|-2 „ 

P\m\ - -7. -^- n -^ . . . = 0 , 

P | m |-|-1 P \ m\+2 


(64) 


or one of its inversions. 


2. Radial solutions 


Let R 2 = \/AR+ and Ri = R-, so that the radial equations flTfll - fl38l) become 

f d iK\ 

[ d iK — MAr\ X — ifir 

U +-A-= 

Near the horizon the ‘ingoing’ solution goes as, 

lim R+ ^ (r — lim R_ —>• (r — r+y, 


(65) 

( 66 ) 


(67) 
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and at infinity the decaying solntion resembles 


where 


lim r'^ lim 

r^CxD r^CxD 



z/ = M 


2MiCjr^ 
r+ — r_ 
— 2u^ 

q 


? 


(? ^ 


( 68 ) 


(69) 

(70) 

(71) 


We may write the solntion in terms of a series, 


R_ 

(r - r+)/ 2 + 


CXD 


= (r — r-^y{r — r_) 


—cr+uqr 




e- 7 iSfc 
fc =0 


r — r+ 




(72) 


Inserting Eq. fl72|) into the radial eqnations leads to a three-term matrix-valued recnrrence 
relation, 


+ /3o^o — 0 , 

^fc^fc+l “ 1 “ f^k^k R 'Yk^k—l O 5 


Here, the matrices are of the form 


ttfel Oik2 
0 0(^4 


ttfc = 

and the matrix coefficients are 


f3k = 


fdkl fdk2 
(dk3 l^kA 


k>0. 


Ik = 


7 fci 0 
7fc3 7fe4 


ttfci = {k a -\- l)(r+ — r_) — i{2Mujr^ — am), 
ak 2 = -i\ + z/ir+)(r+ - r_), 

= (fc -|- CT -|- — r_) -|- i{2Mujr^ — am), 


(73) 

(74) 


(75) 


(76) 

(77) 

(78) 


/Sfci = (r+ — r_) (q'(r+ — r_) — 2{k -\- a) -\- v) -\-2i {2a^u — am) , (79) 

/dfc 2 = (A-hi/nr_)(r+-r_), (80) 

/?fc 3 =-(A-i/ir+)(r+-r_), (81) 

PkA = {r+ — r_) (g(r+ — r_) — 2{k a) u 1) — 2i (2a‘^u — am) , (82) 


Iki = {k a — 1 — z/)(r+ — r_) — i{2Mur_ — am), (83) 

7fc3 = (A - i/ir_)(r+- r_), (84) 

7 fc 4 = (A; -I- (T — I — z/)(r+ — r_) -|- i{2Mojr_ — am). (85) 
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ations can be solved using matrix-valued contin- 
|. We seek the roots of M^q = Oj where 

M = /3o - ao [/3i - Oil {(32 + 0 L 2 M) 72 ]”^ 7i , (86) 

and 

■^n {(3n+l “1“ ^n+l-^n+l) 'Yn+l ; (87) 

with denoting the matrix inverse. Non-trivial solutions exist if 


Matrix-valued three-term recurrence re 
ued fractions, as described in Refs. 58 


det |M| = 0 . 


( 88 ) 


We used a numerical root hnder to hnd pairs that simultaneously satisfy Eq. flM)) 

and Eq. (1881) . As initial estimates for the root-hnding algorithm, {AojWo}, we made use of 
the series expansions in Ref. 78| and the hydrogenic approximation, Eq. (]3|). 


B. Time domain evolution 


Several groups have investigated the excitation of bound states by generic initial data, 
using time-domain codes. The ma iori ty of work has focussed on bosonic helds, in the scalar 
54, 84-^ and Proca cases . The evolution of Dirac bound states was investigated 




53| 


by Zhou et al. in Ref. 

We developed a 1-l-lD time-domain code to solve the coupled radial equations for the 
Schwarzschild case, written in the form. 


at or^ r 


^ _ dF 
dt dr^ r 


F - 


(89) 

(90) 


where F = Ri + R 2 and G = i(i ?2 — Ri) and A = ... , — 2 , — 1 , 1 , 2 ,.... We used the method 
of lines with a fourth-order Runge-Kutta integrator, and fourth-order hnite differencing on 
spatial slices. To implement the absorbing boundary condition at the event horizon, we 
used the ‘Perfectly Matched Layer’ method, which was previously applied in the black hole 
context in Ref. 2^. This entailed adding extra terms —^{r^)F and — 7 (r*)G' to the equations 
above, where 7 (r*) is a non-negative smooth function which is zero for r* > —150M. At 
the far-held boundary we imposed a rehecting boundary condition by setting F(rniax) = 0. 
We moved this boundary to a large radius, rmax ~ 4000M, to suppress its effect on the 
excitation and evolution of low-lying bound states. To analyse the spectrum, we computed 
the ‘spectral power’ 

P(a;) = |F|2 + |g|2 (91) 


from the square magnitudes of the (discrete) Fourier components F{u) and G{u) 
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IV. RESULTS 


A. Perturbation theory: <C 1 


Many authors have shown that, in the limit M/i -C 1, the bound state spectrum is 
hydrogenic; see Eq. (|3]). In this regime, the spectrum is degenerate, as it depends only on 
the principal quantum number n = £ + h + 1 (where ^ is the orbital angular momentum 
[i = •^ + 5', where S G {—1/2, +1/2}] and h = 0,1,... is the excitation number), rather than 
on j, £ and h individually. 

How are degeneracies broken at higher orders in M/i? The spectrum may be written as 
a series in M/x, as follows. 


Ur //i 1 + +- 


(92) 


The ‘hydrogenic’ term E//’^ = 




’( 1 ) 


QiiMfiY 


is familiar; the ‘hne-structure’ term 

and ‘hyperhne-structure’ term E^^^j ~ O ((M/i)®(^)) are less well-known 37, MI, l88|. Here 
we review these terms, and test against numerical data. 


1. Fine structure 


In the case of the hydrogen atom, the hne-structure correction arises at O(a^) where a 


is the hne-structure constant, and it is given by 


89 


£7(1,Hyd) ^ ^ A 


(93) 


The second term, due to a spin-orbit coupling, breaks the degeneracy between states of 
the same orbital angular momentum but opposite spin projection (i.e. the j = £ + 1/2 and 
j = £ — 1/2 states, for £ > 1). 

The hne structure of the Schwarzschild bound state spectrum was calculated in Chap. 5 
of Ref. 88||. The calculation is outlined in Appendix |Dl where a small but vital correction 


is identihed. The key result is 


y{l) _ 

‘^jin 


{MfiY A5 3n 




3n 


2j + l 2i+l 


(94) 


The hnal term, not present in the hydrogen case ([93]), implies that the bound-state spectrum 
depends on three quantum numbers individually: j, i, n (rather than just j, n as in Eq. fl9^ 1. 

Figure [1] confronts Eq. (IMD with numerical data. Here, we plot the diherence between 
ur/^ and the hydrogenic result 1 + En \ and we scale this diherence by The 

plot shows excellent agreement between the prediction of Eq. (j94jl (indicated by horizontal 
dotted lines) and the numerical data (thicker lines), for various modes. This agreement in 


17 


















FIG. 1. Fine structure of Dirac energy levels for Schwarzschild black hole. The solid lines show 
numerical data for the difference Re(a;/^) — (1 — where n = l + n + 1 is the principal 

quantum number, with the y-axis scaled by {M The thick lines show the ground-state 
modes (h = 0) with A = —1 {j = 1/2, ^ = 0), A = -t-1 (j = 1/2, ^ = 1), A = —2 (j = 3/2, 
£ = 1), X = +2 {j = 3/2, £ = 2), A = —3 {j = 5/2, i = 2). The horizontal dotted lines show the 
fine-structure prediction of Eq. ([Ml) , with coefficients —21/8, —25/8, —13/8, —87/40 and —57/40, 
respectively. 


the Mfi —0 regime increases our confidence in the validity of both the numerical method. 


and the perturbation-theory analysis t 
Many years ago, Ternov and Gaina 


na' 


leads to Eq. 


41| also calculated hne-structure splittings for bound 
states, by applying standard perturbation theory techniques. Unfortunately, their hne- 
structure result at 0((M/i)^) is not found to be in agreement with Eq. (1941) . nor with 


our numerical data. We suggest that the calculation in Ref. jdlj is invalid due to the 


coordinate system singularity in Schwarzschild/Boyer-Lindquist coordinates, which causes 
physical quantities to diverge as r ^ r+. We note that, by contrast, Eq. (IMD was derived 
with a horizon-penetrating coordinate system. 


2. Hyperfine structure 

In the hydrogen atom, the interaction between the magnetic dipole moments of the 
nucleus and electron leads to a spin-spin splitting of the £ = 0 state, with the anti-aligned 
spin conhguration lying at a lower energy. This is an example of ‘hyperhne’ splitting, as it 
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FIG. 2. Hyperfine structure of Dirac energy levels for a slowly-rotating Kerr black hole. The plot 
shows numerical data for the difference in the energy of corotating and counter-rotating states {m = 
+j and m = —j), after rescaling by ^{M/a)n^{M. The solid lines show the a = O.OIM dataset, 
and the points show the a = 0.02M dataset; their agreement is evidence that the interaction is 
linear-in-a at this order. The data supports the result of Ref. [371] that rotation-induced hyperfine 
splitting scales with {am/M){Mn)^ at leading order (see text). 


arises at subdominant order a^(me/mnuc.)- The hyperfine transition between £ = 0 levels 
generates a 21cm neutral hydrogen line, of importance in astrophysics. 


In Ref. 4l||, it was calculated that the rotation of the black hole leads to a ‘hyperfine’ 
splitting at O ((Mp)^(^)). It was asserted therein that only the modes i > 1 are split i 
this way. 


m 


In Fig. H] we attempt to extract the hyperfine correction from our numerical data, in 
the slow rotation regime (a < 0.02M). The plot shows the difference between the energies 
of maximally co- and counter-rotating modes, after rescaling by |(M/a)u®(Mp)“^. We 
find that the co-rotating mode is more weakly bound than the counter-rotating mode, as 
expected. The numerical data suggests that the scaling of Ternov et al. |4]| is indeed correct. 
However, we hnd that the £ = 0 mode is also split at this order, and the numerical coefficients 
found in Eq. (38) of Ref. [^ are not consistent with our data. 


19 













3. Decay rates 


In Ref. 37j , Ternov et al. used an asymptotic matching method to derive a key result for 
the imaginary part of the frequency of Dirac bound states on the Kerr-Newman spacetime 
in the regime Mfi/j -C 1, am/M -C 1 (see Eq. (29) in Ref. 37|]). In the Kerr case, 


(U/Z/i -a jin {M^) 


i+2l+2S 


(96) 


where S = ±1/2 and 


^jln 


r+ 


r_ 


± r_ 

i+i/2 r 


1+21 


r_ 


X 


n 

p=i 


i± 


(2/ ± l)(r+ ± r_) 
4r2 


2S 


(n±£)! 


^4+2£(2£)!(2£ + i)!(n-£- 1)! 


(p-1/2)2 


(96) 


with r = (2Mr+p — am) / (r_|_ — r_). 

Figure [3] shows numerical data for the imaginary part of bound states of the Schwarzschild 
black hole, for n = 1... 3 modes. The plot illustrates how the decay rate has as a power-law 
scaling in the M/i -C j regime. The leading-order results of Eq. (l95|l are shown as dotted 
lines. These results are found to be in good agreement with the data in the small-M/i 
regime. 


B. Time-domain evolution 


Are bound states typica 


number of studies 



excited by generic initial data? For the scalar held case, a 

iJ 


86 


have concluded that the answer is affirmative (see 
Ref. ^ for a Green’s function analysis of the excitation factors in the initial-value formu¬ 
lation). Thou gh the conclusion is expected to be similar in the Dirac case, it has received 
less attention and we briehy investigate it here. 

Figure Hlshows the power spectrum P{u:) [Eq. fl?I|) ] resulting from a typical time-domain 
evolution, with parameters Mfi = 0.2, a = 0, A = —1, and Gaussian initial data F{r) = 0, 
G(r) = exp (r* — r®)/(2cr^)j, with = 90 and a = 40. There is clear evidence here 
that the hrst three bound states were excited by the initial data, at the expected frequencies 
{ojk/h ~ 0.9741726, 0.9940106, and 0.9974677). Furthermore, the ‘spectral lines’ have the 
expected Lorentzian prohles, P ~ l/[(i^ — i^r)^ ± i^/]- The widths are found to be in 
proportion to the imaginary parts of frequencies {—ui/fi ~ 3.755 x 10“^, 4.979 x 10“^ and 
1.407 X 10 ®). Gonsequently, the spectral line for the ground state is found to be the widest, 
as this mode decays most rapidly. 
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FIG. 3. Power-law scaling of the decay rate of bound states on Schwarzschild spacetime. The 
solid lines show numerical data for —Im(a;///) for the first two modes with A = —1 (j = 1/2, t = 0), 
A = -|-1 (j = 1/2, £ = 1) and A = —2 (j = 3/2, £ = 1). The dotted lines show the asymptotic 
results and (see Eq. (I9^ h 

C. M/i ^ 1 regime 

In the ‘semi-classical’ regime M/i S> 1, in which the gravitational length scale is much 
longer than the quantum-field length scale, we expect the key features of the spectrum 
to relate closely to the properties of timelike geodesics. Indeed, this quantum-to-classical 
correspondence emerges via a standard WKB analysis, as shown in Appendix lAl 

For circular geodesics on the Schwarzschild spacetime, the (dimensionless) energy Cj is 
related to the angular momentum L and orbital radius fo by: 



The orbital radius fo is given in terms of L in Eq. (IA7p . For L ^ 1, a) ~ 1 — \L~‘^ — -|- 

.... In Appendix 13 it is shown that the geodesic parameter L should be associated with 
{£ + l/2)/M/i for a scalar field. 

Plot (a) of Fig. [5] shows numerical data for the h = 0 energy level of a selection of angular 
modes in the range —4 > A > —10. The energy level is shown as a function of M/i/|A|. We 
may associate this parameter with l/L in the classical limit. Making this association, the 
geodesic energy level of Eq. fl97l) is shown as a solid black line. The agreement between the 
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FIG. 4. Power spectrum for a time-domain evolution of a massive Dirac field {Mfj, = 0.2) on 
Schwarzschild spacetime (a = 0). The j = = 0 (A = —1) mode was evolved with the 

1-l-lD scheme outlined in Sec. IIVBI starting with generic initial data described in Sec. IIVBI up to 
imax = 4 X lO^M. The plot compares numerical data for the power [Eq. (I9ip . red symbols] with the 
spectral lines for the first three bound states [black vertical lines] determined via the frequency- 
domain method of Sec. lIII Al) . The expected Lorentzian half-widths are shown as horizontal arrows. 


data and geodesic prediction is evident. A dotted vertical line indicates the innermost stable 
circnlar orbit (rg = 6M, L = l/\/T2), beyond which stable circnlar orbits do not exist. Our 
numerical data suggests that, though bound-state solutions exist beyond this limit, they are 
rapidly-decaying. 

Plots (b) and (c) of Fig. Oshow the imaginary part of the frequency of these modes, as 
a function of M/r/|A| 1/L. Plot (b) shows that the decay rates are very similar at the 

position of the innermost stable circular orbit, L = \/T2 (the point at which the angular 
momentum potential barrier separating the circular orbit from the horizon disappears). Plot 
(c) shows that the decay rate is consistent with 


jj) ~ exp (^I3{L)M 


(98) 


where /5(h) is some quantum-tunnelling factor, determined from an integral across the po¬ 
tential barrier. 
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FIG. 5. Bound state frequencies in semi-classical regime {M^. S> 1, |A| 1), for a Schwarzschild 

black hole. The dotted line at M^/\\\ = (12)“^/^ corresponds to the innermost stable circular 
orbit at r = 6M. 


D. Kerr bound states: decay 

Now we focus attention on the decay rate of modes on a rotating black hole spacetime. 
Figure [6] shows the imaginary part of bound state frequency as a function of M/i, for both 
co-rotating (m > 0) and counter-rotating (m < 0) cases. For low couplings M/i -C 1, we see 
that the counter-rotating case is rather similar to the Schwarzschild case: the decay rate is 
governed by a power law, in accord with Eq. ([95]). 
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M|i 


M|i 


j=3/2, 1=1 



M|i 


FIG. 6. The decay rate of bound states on Kerr spacetime for co-rotating {m = +j) and counter¬ 
rotating {m = —j) modes. The plots show —Im(t(j//i) for the j = 1/2, £ = 0 [left] j = 1/2, 1 = 1 
[right] and j = 3/2, £ = 1 [lower] states, as a function of M/r, on logarithmic axes. 

The co-rotating modes exhibit an interesting feature. In the regime u < mfln the decay 
rate is suppressed; whereas outside this regime, the decay rate increases rapidly. The plots 
suggest that, for rapidly-rotating black holes a > 0.99, the decay rate has a local minimum 
precisely at ca = mflH. 

Figure [7] shows the decay rate as a function of cu/mflH, where Qh is the angular frequency 
of the event horizon. The plot makes it evident that the local minimum coincides with 
the ‘critical’ superradiant frequency ujc = mQH- The data suggests that minimum value 
decreases with increasing a, which hints at the possibility that there may be a non-decaying 
‘critical’ mode in the extremal case a —)■ M. This possibility remains to be investigated. 


E. Kerr bound states: spectrum 

Now we turn attention to results for the Kerr energy spectrum. In the limit M/i ^ 0 the 
spectrum is well-understood, as outlined in Sec. IIV Al but outside this regime the spectrum 
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0.01 


j = 1/2, l = 0 
j = 1/2, 1 = 1 



FIG. 7. Decay rate for co-rotating modes m = +j at o = 0.999M as a function of ojIuJc.-, where 
ujc = ttiVLh is the critical superradiant frequency. The plot shows that decay is suppressed for 
u)u) < 0 , and that a local minimum arises at cu = Uc- 


has not been examined in detail. 

Figure [8] shows the spectrum of the lowest modes n = 1 and n = 2 (where n is the 
principal quantum number) for a Kerr black hole with a = 0.9M, as a function of Mfi. Let 
us make several observations. (1) At low M/i, the j = 1/2 co-rotating (m = -1-1/2) state is 
higher-lying, and decays more slowly, than the counter-rotating state (m = — 1 / 2 ). ( 2 ) As 
Mfi increases, the roles are reversed, with the co-rotating state lower-lying for Mfi > 0.2 
and faster-decaying for Mjj, > 0.33. (3) At higher couplings, the j = 1/2 states appear 
to achieve large binding energies; but here the states are highly transient (with no classical 
analogue) and may not be physically signihcant. (4) At low M/i, the ^ = 1 states are ordered 
according to the hne and hyperhne structure equations: the j = £ — 1/2 states [solid] are 
lower-lying and faster-decaying than the j = £-1-1/2 states [dotted], and the counter-rotating 
states (m < 0 ) are lower-lying than the co-rotating states (m > 0 ). Fig. [Hl[b) shows how 
this hierarchy changes as M/i increases, and the states become more transient. (5) The 
maximally-corotating modes (e.g. m = j = 3/2, £ = 1) are slowly-decaying in the regime 
Re(ci;) < mfln- 

In Table H we present a sample of numerical data for bound state frequencies, together 
with the corresponding angular separation constants. 
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FIG. 8. Spectrum of n = 1 and re = 2 Dirac bound states for Kerr black hole at a = 0.9M. 
Plot (a) shows the real part of the frequency, Re(a;)//r, as a function of M/r, for j = 1/2, ^ = 0 
(dashed), j = 1/2, £ = 1 (solid) and j = 3/2, i = 1 (dotted) modes. Plot (b) shows the fine 
structure, [Re(a;)//x — l]/(M/r)^, for the re = 2 states. At low M/r, the j = 1/2, £ = 1 states [solid] 
are lower-lying than the j = 3/2, £ = 1 states [dotted]; and the counter-rotating states (rre < 0) 
are lower-lying than the co-rotating states (rre > 0). The ordering changes as M/r increases, as the 
imaginary part of the frequency increases: see plot (c). 

F. Kerr bound states: wavefunctions 

The components of the Dirac current in the ‘ingoing Kerr’ coordinate system were con¬ 
sidered in Sec. Ill B 81 In obtaining (l5Tll - (]55|l , we tacitly assumed that the frequency uj is 
real. For bound states as defined in Sec. HU for which this is not the case, Eq. (l55ll should 
be multiplied by an additional factor of exp(2Im(a;)t). After integrating over the two-surface 
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Mode m 

Re(w)//x Im(a;)//x 

Re(A) Im(A) 

7 = 1/2, ^ = 0, -1/2 
+1/2 

0.95120580 -3.1293497 x lO"^ 

0.92200086 -2.3193867 x lO-^ 

-1.2612238 +5.6180382 x 10"® 

-0.7440728 -4.1963431 x 10"® 

j = l/2,i = l, -1/2 
+1/2 

0.98580062 -8.1786335 x lO"'^ 

0.98662631 -6.3064813 x lO"*^ 

+1.1074300 -1.6229636 x 10"^ 
+0.9348739 +9.8118563 x 10"® 

7 = 3/2, £ = 1, -3/2 
-1/2 
+1/2 
+3/2 

0.98744302 -3.2293333 x lO''^ 

0.98767652 -4.1491804 x 10"® 

0.98785317 -7.7639988 x 10"® 

0.98798850 -4.2000702 x 10"® 

-2.2672880 +6.9736117 x 10"® 
-2.1091337 +3.7495425 x lO”’^ 

-1.9333717 -3.6997312 x lO-^ 

-1.7325948 -9.0746837 x lO-^ 


TABLE I. Sample numerical data for Dirac bound states (n = 1,2) for a Kerr black hole at 
M/x = 0.3, a = 0.9M. The columns show the real and imaginary parts of the frequency, co/fj,, and 
the angular separation constant, A. 


of constant r, and applying the normalization condition flT2ll . we have 


(T(o) + T(i)) 


where 


T(0) = 


+ 2Mr + a^) | 


(99) 


( 100 ) 


[r — r_j 

with T(i) associated with the second term in Eq. ([55]). These terms are hnite at the onter 
horizon. 


In the limit M/x —)■ 0, the radial solntions are hydrogenic in form ^ (see Appendix |D|) . 
The radial fnnctions may be written in terms of Lagnerre polynomials in the dimensionless 
variable x = [MpY {r/M)rr^. The radial fnnction has n nodes, where fi = n — ^ — 1 is the 
excitation nnmber. 

Fignreinishows the radial prohle of T(o) for the j = m = 1/2, i = 0 mode at M/x = 0.38, 
for a variety of Kerr parameters a. Changing a modihes the shape of the prohle somewhat. 
In this case, the peaks in the prohle move closer to the horizon, and the valne on the horizon 
also increases, as shown in the inset. 


V. DISCUSSION AND CONCLUSIONS 

In this article we have explored the spectrnm of bonnd states of the massive Dirac eqna- 
tion on Kerr spacetime: ‘trapped’ modes which are ingoing at the horizon and fall oh 
exponentially towards spatial inhnity. Let ns now briehy review the key resnlts. 

In Sec. ini we formnlated the Dirac eqnation on Kerr spacetime, reviewing (i) the sep¬ 
aration of variables, (ii) the absence of snperradiance, and (hi) the violation of the weak 


27 













6 



FIG. 9. Profile of probability density from Dirac current, for j = 1/2, n = 3 bound states on 
Kerr black holes with M/x = 0.38 and a/M = 0, 0.5, 0.9 and 0.99. The plot shows T^g), defined in 
Eq. (jlOOl) . which is related to the Dirac probability density in the ingoing Kerr coordinate system. 
The inset shows that the density is finite at the outer event horizon, as expected. 


pri nciple. Eschewing the Newman-Penrose 2-spinor formalism used in p ioneering 

with 


63 


64, 1 ^ , we favoured instead the 4-spinor formalism 0 , 0 , 0000 , 


energy 
works 

the Weyl representation; positive spacetime signature; and Carter’s canonical (‘symmetric’) 
tetrad. 

In Sec. Uni we presented a practical method for computing the bound state spectrum, 
reducing the problem to that of Ending the roots of certain (matrix-valued) three-term 
recurrence relations. We also briefly described a time-domain method for the Schwarzschild 
case, with scope for extension to the Kerr case. 

In Sec. IIV Al we confronted small-Mp asymptotic results, obtained in the 1980s by Ternov, 
Gaina and coworkers with new numerical data. We concluded that (i) the decay of 
bound states is well described by Eq. fl95l) : and (ii) the scaling of the fine-structure and 
hyperfine-structure terms has been correctly deduced. However, the numerical coefficients 
for the (hyper)fine-structure terms were found to be inconsistent with the data. We found 
that, instead, our data is fully consistent with the fine-structure result presented in Eq. 


and Appendix [Dl 


which was derived in Ref. 

In Sec. IIVBI we used a time-domain code to demonstrate that bound states are typically 
excited by generic initial data. In Sec. IIV Cl we briefly considered the large-M/i regime, 
showing in Fig. [5] that our new results are consistent with ‘semi-classical’ expectations from 
a WKB analysis (e.g. Appendix |A]) . In Secs. IIVDfHVFl we examined the novel features of 
bound states of rapidly-rotating black holes. We found that the decay rate of the maximally- 
corotating mode is strongly suppressed in the regime oj < mfln (Figs. ISHZl)- 




























Let us now discuss some of the implications of our findings. For the scalar held, modes 



stationary, due to the lack of superradiance. On the other hand, as decay may be strongly 
suppressed (Fig. [7]), it is plausible that rapidly-rotating black holes can support very long- 
lived co-rotating Dirac hair. 

Let us now highlight several avenues for future work which could lead to a more complete 
understanding. Despite the progress described above, this work lacks: (i) correct coefficients 
for the hyperhne structure at 0{{am/M){M (ii) an asymptotic result for the rate 
of decay for co-rotating modes for uj < mfln, a < M; (iii) asymptotic WKB results at 


subleading order (for M/i S> 1) to account for the effects o 
analysis of the Dirac equation in second-order form 




spin. Revisiting the asymptotic 
could pay dividends. 

Another possibility would be to examine the bound-state wavefunction inside the black 
hole horizon. It is well-known that the inner (Cauchy) horizon at r = r_ has the effect of 
repelling geodesics. Therefore, it seems unlikely that the bound states investigated here will 
be ‘ingoing’ at r = r_, in general. Nevertheless, there may exist particular frequencies for 
which the mode may be ingoing at both horizons. This remains to be investigated. 

A further topic for investigation is the excitation of Dirac bound states by generic initial 
data on Kerr spacetime. Two approaches suggest themselves. First, an extension of the 
Green’s function analysis of Ref. ^ to the Dirac case. Second, an extension of the time- 
domain approach of Sec. IIIIBl and Refs. [2J, to the a > 0 case. 

In this work we have treated the Dirac held as a ‘classical’ held, which is of course 
not the case. Unruh’s second-quantized analysis 79|] (1974) showed that fermions, as well 
as bosons, experience an instability to spontaneous particle creation (Unruh-Starobinski 


radiation [79|, [9^); that is, a quantum version of superradiance. In 1983, Gal’tsov et al. 
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considered the hlling of Dirac bound states in the Schwarzschild case, concluding that the 
distribution is thermal {N ^ [8(1 -|- exp{u}/ kTh))]~^) in the limit M/i -C 1, where Th is 
the Hawking temperature. The hlling of states on the Schwarzschild spacetime was also 
considered in Ref. 46(]. Hartman et al. 47|] have considered the Kerr case, arguing that 
all bound states with energies u < mVtu will be hlled in Unruh’s vacuum, creating a 
(stable) ‘Kerr-Fermi sea’ which extends outside the ergosphere. This intriguing possibility 
undoubtedly deserves some further consideration. 
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Appendix A: WKB analysis of Schwarzschild bound states 


Consider a scalar field $ on Schwarzschild spacetime, satisfying □$ — = 0, which 

may be decomposed into modes in the standard way, 

1 


$ = -e—*n(r)Fz™(0,0). 


The radial eqnation is 




— Ten - EoR) = 0, 


/ \ it \ I 2 I T 1) 2iff 
Voir) = fir) . + — 


(Al) 


(A2) 


where /(r) = 1 —2M/r and r* = r + 2Mln(r/2M —1). In the regime Mjj, S> 1 (i.e. S> Ac), 
this may be written in dimensionless form as 

^ '1 

+ (M/i)2 C(f) U = 0, f/(f) = - ■P(f) + O ((M/i)-2) , (A3) 




and 


V[f) = fir) 1 + ^ 


(A4) 


where f = r/M, f* = r^/M, Cj = a;//i and !/ = (£ + l/2)/(M/i). Note that (IA3P resembles 
a Schrodinger eqnation with a large parameter M/i \/2m/h, and thns we may apply 
standard WKB methods to reach the Bohr-Sommerfeld condition for bonnd states: 
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f/V2^A = ^ (n + 1/2), n e N. 
M/i 


(AS) 


At this point we note the close analogne between [/(f) in flA3p and the right-hand side 
of the geodesic ‘energy’ eqnation. 


= [/geo(r) ^ 


/(f) (1 + ^ 


(A6) 


where 8 = ft and C = r‘^(j)/M are constants of motion, and the overdot notation indicates 
differentiation with respect to proper time. Circnlar orbits are defined by dU^eo/dr = 0 with 
stable orbits satisfying d^f/geo/dr^ < 0. The former condition leads to an eqnation for the 
circnlar orbit radins in terms of the dimensionless angnlar momentnm: 


fo = W[l+Vl-l2/L 


(A7) 


The latter condition implies that fo>6 and > 12 for stable orbits. After expanding 
aronnd the circnlar-orbit radins, and inserting into flASp . we obtain 


0,2 = + ..., 


(AS) 
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where Vq = V(fo) and V" = 


d?y 

~dr^ 


. Equivalently, 


ro 


^ t.1/2 , IV" ,, ,{n + l/2) , 

M/r +••• 


(A9) 


In the weak-held regime, where fo ~ 3 > 1, Vq ~ 1 — fg V" ~ 2/f^ we reach Eq. 

(after reinserting dimensionfnl constants). 

Appendix B: Carter tetrad and spin connection 


The inverse components of Carter’s tetrad are 


e^„dx^ [dt 


asm 


eldx^ 


P 

sin 9 
P 


[—adt + (r^ + a^)d(j)^ , 


e\dx^ =-^dr, 
\/A 

eidx^ =pd9, 


(Bl) 


Noting the anti-symmetry relation ujfj_ab = the non-trivial components of the spin 

connection are 




M(r^ — cos^ 9) 


2aMr cos 9 


P^ 


^r03 — 


ar sin 9 
p'^'/a ’ 


^eo3 — 


i\/A cos 9 


P 

a sin^ 9 
= — Za —B, 


<^(1)13 


r-v/AsinP 

9 ’ 


(^t23 — 


^rl2 — 


(^012 — 


(^if>02 


1^023 


P 

sin 9 cos 9 


,2 

r\/A 


P^x/A ’ 


P^ 


i\/A sin 9 cos 9 


P" 


cos 9 
P^ 


A, 


where A = p^A + 2Mr{A -|- a?) and B = (Ar cos^ 9 — AM cos^ P -|- -|- Mr"^. 

The matrices T^, dehned in Eq. flT^ . are given by 

M (Q~‘^ (Ts 


T, = 


T. = - 


2 \ O 
1 / asinP 


\/A 


r. ^ -1 (VA) 


= -\/AsinP 


O 

-0*~^ (^3. 


Q ^2 

o 

(ip)~^CT2 

O 

{ig)~^ ai 


O 

-g*~^a2 


O 


■{ig)* (72 

O 

o -(ip)*”Vi 

where zu = i cos 9 — ag~'^(g -|- M) sin^ 9. 


+ 


1 (zua^ O 
O —w* (J 3 / 


(B2) 


(B3) 
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Appendix C: Connection with Finster et al. 


After multiplying Eq. (I5S]) by a^, we may write once more in the four-matrix form 


(7^(r) + Qie)) 


v~ 


= 0 


(Cl) 


where 

/ i/ir 

0 

\/AD+ 0 \ 


( a/r cos 6 

0 

0 


n{r) = 

0 

\/AD_ 

—ijir 

0 

0 \/AT>_ 

—i/rr 0 

, Q{r) = 

0 

0 

—a/i cos^ 

a/i cos 9 

0 

0 


V 0 

\/ad+ 

0 ijir y 


V 

0 

0 

—a/i cos 









(C2) 

with 












V± = ± ^ -|- a^) -|- am] 



(C3) 




C± = de + ^ cot 6 ± {m esc 6 — 

aoj sin 6). 



(C4) 


A spinor transformation takes this to the form found in Refs. 
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Appendix D: Fine-strncture calculation for Schwarzschild bound states 


Here we give an overview of a calculation in Chap. 5 in Ref. 88|, which leads to the hne- 
structure result, Eq. (|M|l . The calculation starts with the Dirac equation in the ‘Newtonian’ 
gauge (PainlevAGullstrand coordinates) which can be written in Hamiltonian form 
idt^J = Hip where 


Hip =+ fij^ip + Hiip, Hjip = Ip, (Dl) 


with gamma matrices from the Dirac-Pauli representation. This resembles a flat-space equa¬ 
tion with a novel interaction term. Applying a Foldy-Wouthuysen transformation 97| leads 
to 

H = j^fi + Ho + Hi + (D2) 


where 


Ho = —i^O‘^ + S, H,= 


1 .. 1 
- 




0 ,£ 


(D3) 


2/i' ' 8/^3' 8/i2 

and the Pauli-even and odd operators are O = and S = Hj. The Foldy-Wouthuysen 

transformation renders the Dirac equation in block-diagonal form, with upper (lower) com¬ 
ponents representing particles (anti-particles). 
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The zeroth-order equation idt'ip = + Hq'i/j can be written in more familiar form by 

introducing the phase-transformed wavefunction, 

leading to 


{idt - /i)T = — 

2/i r 


Introducing a separation of variables, 




(D5) 


(D6) 


where x is any constant two-spinor, leads to a time-independent ID Schrodinger equation 
with a 1/r ‘Newtonian’ potential. 


1 f 2d + 




R 


'Tli 


Mjj, 


R. 


■nl' 


(D7) 


2/i \dr‘^ r dr 

This equation has standard hydrogenic solutions for the bound state energy levels 8n^ = 
— (M/i)^/2n^ and wavefunctions, 

Rni{r) = (D8) 

where x = (2/?i)(M/i)^(r/M), is an associated Laguerre polynomial, Yi^ is a spherical 
harmonic, and Ane is a normalization constant. 

To deduce the hne-structure correction to the energy, one calculates the expectation 
values of terms in the first-order Hamiltonian Hi, Eq. (ID3H . when closed with the zeroth 
order wavefunctions, that is. 


A,] = 




,,iMfiy/Sr/M-qj 


nl 


The details of this calculation may be found on plOO-105 of Ref. 


correction is made up of two parts, 

{M^^y 


7°C»^ ) = 


8/i^ 


8n‘^ 


-15 + 


48?7, 


+ 


9n 


(D9) 

The fine-structure 


(DIO) 


and 


8/i3 


O, 0,8 


Here we have corrected an error in Ref. 


(this error may be traced back to the final entry of Table 5.1 of Ref. 
be a ‘2’). Taking the sum of contributions leads to Eq. 


2£ + l ' £(£+l)(2£+l) 

(M,.)- n(12A + 3) 

8ii* lit + l)(2t + 1)' ' ’ 

in the coefficient of the second term of Eq. flDlOp 

, where the ‘4’ should 
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